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O ' ABSTRACT 
O 

We use a complete sample of active galactic nuclei (AGN) selected on the basis of 
relativistically beamed 15 GHz radio flux density (MOJAVE: Monitoring of Jets in AGN 
with VLBA Experiments) to derive the parent radio luminosity function (RLF) of bright 
radio-selected blazar cores. We use a maximum likelihood method to fit a beamed RLF 
to the observed data and thereby recover the parameters of the intrinsic (unbeamed) 
RLF. We analyze two subsamples of the MOJAVE sample: the first contains only 
^ . objects of known FR II class, with a total of 103 sources, and the second subsample 

adds 24 objects of uncertain FR class for a total of 127 sources. Both subsamples 
exclude four known FR I radio galaxies and two gigahertz-peaked spectrum sources. 
We obtain good fits to both subsamples using a single power law intrinsic RLF and 
a pure density evolution function of the form z m exp [— l/2((z — zo)/a) 2 ~\ . We find 
that a previously reported break in the observed MOJAVE RLF actually arises from 
using incomplete bins (because of the luminosity cutoff) across a steep and strongly 
evolving RLF, and does not reflect a break in the intrinsic RLF. The derived space 
density of the parent population of the FR II sources from the MOJAVE sample (with 
£i5GHz > 1.3 x 10 25 WHz _1 ) is approximately 1.6 x 10 3 Gpc -3 . 



Subject headings: galaxies : luminosity function, mass function — galaxies : evolution 
- galaxies : active — quasars : general — BL Lacertae objects : general - 
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Introduction 



The radio luminosity function (RLF) of active galactic nuclei (AGN) and its redshift de- 
pendence are important quantities in understanding the physics of AGN and their cosmological 
evolution. In the case of AGN selected on the basis of relativistic emission (i.e., blazars), it can 
also provide information about the parent population from which an observed sample is drawn. A 
parameterized luminosity function (LF) can also be useful for producing M onte Carlo simulations o f 
populations to compare with statistical properties of observed AGN (e.g.. iLister Marscherlll997l ) 
as well as to study those properties of AGN that are difficult to observe directly. The intrinsic RLF 
can also be useful for predicting the numb er of 7-ray blazars to be observed by future surveys (e.g., 
GLAST; also see ILister Marscherlll999l ) as well as for determining how rare individual blazars 
are in the general AGN population. 



According to contemporary AGN unification schemes (see review bv lUrrv Padovanilll995l ). 
various observed classes of AGN (e.g., radio galaxies, quasars, and BL Lacs) can be the result of 
different orientations of essentially the same type of object. One can test unification schemes using 
statistical approaches. For example, if BL Lac objects are highly beamed versions of lower power 
radio galaxies, then the number of BL Lacs should be much smaller then the number of p arent radio 



galaxi e s, because BL Lac s are o riente d at a small angle to the line of sight. Previously, lUrrv et al 



(|199ll ) , iPadovani &: Urrvl ( 19921 ) , and lUrrv fc Padovanil (|1995l ) applied relativistic beaming correc- 
tions (e.g., Cohen et al. 20071 ) to the RLF of high power radio galaxies and f ound it to be compatible 
with the observed RLF of a sample of flat-spectrum, radio-loud quasars. iJackson k, Walll (|1999l ) 
proposed a dual-population unified scheme in which (a) the high-power FR II radio galaxies are 
the parents of all radio quasars and some BL Lac- type objects, and (b) moderate-power FR I radio 
galaxies are the parents of the remaining BL Lac-type objects. They tested this model by beaming 
(using Monte Carlo jet populations with a single bulk Lorentz factor) the low-frequency radio data 
and comparing them with high-frequency radio data. 



The MOJAVE AGN sample ([Lister & Homanll2005l ) is the first large, radio-select ed AGN sam- 
ple fo r which jet kinematic and apparent superluminal speed information are available ([Kellermann et al 
20041 : Lister et al. 2008, in preparation). It is complete with respect to relativistically beamed jet 
emission, and therefore provides a unique opportunity to learn about the intrinsic (parent) RLF of 
blazars. The determination of the intrinsic (non-beamed) RLF is complicated, however, by relativis- 
tic beaming and selection effects. The radio emission from an AGN is highly enhanced by Doppler 
boosting if its jet is relativistic and aligned close to the line of sight. A flux density-limited sample 
of AGN will therefore contain not only sources with high intrinsic luminosity, but also sources with 
lower intrinsic luminosities whose flux densities are Doppler boosted because of their orientation. 
The effect of Doppler beaming on the observed RLF was first calculated for sing le Lorentz factors 
( Urry &; Shaferl I1984T ) and later extended for distributions of Lorentz factors ([Urrv Padovani 



1991 



). ILister! (|20o"j ) extended these studies by deriving fully analytical expressions for the Doppler 



factor distributions and beamed RLFs. 
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Previous studies (see, e.g.. |Padovani fc Urrvlll992l ) started with the assumption that the intrin- 



sic parent LF was that of the FR II galaxies, and then applied beaming in order to compare it with 
the LFs of the steep spectrum radio quasars (SSRQ) and flat spectrum radio quasars (FSRQ). Since 
their sample was selected at a much lower radio frequency (2 GHz), it was subject to contamination 
by sources with large extended emission, as well as uncertainties in the classification of the sources 
into SSRQ and FSQR based on a somewhat arbitrary spectral index cutoff of a r ad = —0.5. In addi- 
tion, in order to simulate the b eamed LFs of SSRQ and FSRQ it was necessary to know the ratio of 



the core to extended emission. iPadovani fc Urryl (jl992l ) assumed a linear relationship between the 



beamed and unbeamed luminosities, but found that they needed different factors to produce good 



fits to the SSRQ and FSRQ. A somewhat similar approach was used by iJackson Sz Walll ([1999) 
in determining the beaming models of the parent populations in their dual-population unification 
scheme. 

The MOJAVE sample is different from previous samples in that it is selected on the basis of 
(highly variable) radio flux densities at a high frequency of 15 GHz, thus effectively eliminating 
contamination from extended source emission. The uncertainties surrounding the spectral index 
cutoff and core-to-extended emission ratios are likewise alleviated. 

In this paper we use the maximum likelihood method to fit a beamed RLF to the observed 
data, from which we recover the RLF parameters of the parent population of the MOJAVE sample. 
These parameters will be used in upcoming studies of the effects of beaming on the blazar properties 
derived from flux limited samples. 

The outline of the paper is as follows: In § [2] we describe the observational sample and our 
method for dealing with incomplete redshift information. In § !3.1l we describe our parameterization 
of the RLF, and in § 13.21 we describe the method used to find the optimized model parameters 
and constraints on the fits. We present the results of the model fitting in § [4] and summarize our 
findings in § [5j 

Throughout this paper we assume (unless stated otherwise) a cosmology with VL m = 0.3, 
VL\ = 0.7, O r = and Hq = 70 km s _1 Mpc -1 . All luminosities are quoted as monochromatic 
luminosities at specific frequency v. We also adopt the following convention for the spectral index, 
a rad : S v oc v a ™ d . 
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Observational Data 



2.1. Sample Description 



The MOJAVE AGN sample ([Lister k Homanl 120051 ) consists of all 133 known bright AGN 
with galactic latitude \b\ > 2.5°, J2000.0 declination greater than —20°, and compact (VLBA0) 
flux density exceeding 1.5 Jy at 15 GHz (2 Jy for sources with 5 < 0) at any epoch between Jan. 
1, 1994 - Dec. 31, 2003. The sky area covered by the MOJAVE survey is 6.00912 sr for the 
northern sky and 2.08012 sr for the southern sky. The sample is selected on the basis of beamed jet 
emission only. The contribution from the large-scale radio emission is effectively excluded by using 
the milliarcsecond scale (VLB A) 15 GHz flux density, since the former tends to be diffuse and has 
a steep radio spectrum. 

Many of the MOJAVE sources exhibit high flux density variability and have been selected based 
on their largest flux density values, th us potentially creating a source of selection bias. However, 
in a detailed study of AGN variability, iListerl (|200ll ) concluded that the effect of the variability on 
the sample selection is small in moderately sized samples because the majority of highly beamed 
sources in the pa rent population (which are preferentially selected in beamed emission-selected 
samples; see, e.g., IVermeulen fc Cohen Ill994l ) would lie well above the survey flux limit and will 
be selected regardless of their flaring levels. Those sources just above the survey limit would be 
statistically balanced by other sources lying just below. 

We present basic data for the MOJAVE sample, including Fanaroff- Riley (FR) and optical clas- 
sificat ions for the MOJAVE sample in Tabled! For optical class we used the lVeron-Cetty &: Veron 
catalog as well as the NASA/IPAC Extragalactic Database (NED). We note that the optical 
classification of bla zars into BL Lacs and optically violen tly variable (OVV) quasars remains con- 
troversial (see, e.g.. lAntonucci Ill993l : iKovalev et al.ll2005l ) but this is not essential to our analysis, 
which is based on Fanaroff-Riley classes. We base our Fanaroff-Riley classification of the MOJAVE 
sources on their extended radio morphology. We assign an uncertain FR class to BL Lacs and inter- 
mediate BL Lac/HPQ sources for which morphological FR classification is difficult (i.e., core, halo, 
unusual morphologies or otherwise lack of prominent hotspots) . For quasars only, in situations were 
morphological classification is difficult we rely on the source's luminosity to assign it a FR class. 
At 178 MHz the Fanaroff-Riley divide occurs at about Li7smh^ ~ 10^ 5,3 WHz" 1 , but this tr ansi- 
tion depends on the host galaxy magnitude (see lOwen Ledlowlll994l : iLedlow &: Owenlll996l ) and 
sources close to this luminosity can b e of either class. The line separating the two classes at 1.4 GHz 
in Fig. 1 of ILedlow Owenl (119961 ) can be approximated as logLpR break ~ 10.2 — |M^ ost over 
the range —25 < M^ ost < —21 (assuming their cosmology of qo = and Hq = 75 km s _1 Mpc -1 ). 
We used this equation to find the FR break luminosity, Lfr break, for a given host magnitude. 
We considered any quasar with an inconclusive morphology to be of FR II class if its luminosity 



The Very Long Baseline Array is operated by the National Radio Astronomy Observatory, which is a facility of 
the National Science Foundation operated under cooperative agreement by Associated Universities, Inc. 
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at 1.4 GHz is at least an order of magnitude larger than LpR break- Unfort unately, host galaxy 
magnitudes are available only for six of the 91 of the quasars in our sample. iPagani et al.l ((2003) 
found that the average absolute magnitude of the host galaxies radio loud quasars at low redshifts 
(z < 0.5) is ( Mr) = —24.0 ± . 5, qp = and Hq = 50 km s _1 Mpc -1 which when converted to the 



cosmology of iLedlow k, Owenl (|1996l ) becomes (Mr) ph —23.1. We assume Mr = —23.1 for the 85 
quasars for which we could not find information in the literature on the absolute magnitudes of 
their host galaxies. We have computed the luminosities at 1.4 GHz from the fluxes found on NED 
(we used the largest listed flux) assuming a spectral index a ra d — —0.7. All of the quasars with 
uncertain morphological FR classifications ended up as FR II class according to this model. Of 
the eight galaxies in the MOJAVE sample, three galaxies (0007+106: III Zw 2, 0415+379: 3C 111 
and 1957+405: Cygnus A) show FR II morphologies, another four galaxies (0238-084: NGC 1052, 
0316+413: 3C 84, 0430+052: 3C 120 and 1228+126: M87) show FR I morphologies, and the galaxy 
2021+614: OW+637 is a gigahertz-peaked spectrum (GPS) source. The other GPS source in the 
MOJAVE sample is the quasar 0742+103. Redshifts are available from NED for all but 12 sources 
(6 optically featureless BL Lacs, and 6 sources without optical identifications). 



According to a contemporary unification scheme (e.g., lUrrv Padovanil Il995l ). the parent 
population of BL Lacs is identified with FR I type galaxies. However, the issue of parent populations 
for BL Lac objects remains under debate. Recent studies o f the host galaxy and extended radio 
emission of radio-selected, l ow energy peaked BL Lacs (e.g., ICassaro et al.lll999l ; iRector fc Stocke 



2001 



Kotilainen et al.ll20 05) appear to rule out the FR I - BL Lac unification scheme in its simplest 



form. 



Jackson &; Walll (|1999l ) proposed a dual-population scheme in which FR II radio galaxies are 
the misaligned parents of flat-spectrum quasars and some BL Lacs. We adopt this unification 
model and exclude 4 FR I galaxies (0238-084, 0316+413, 0430+052, and 1228+126) from the 
sample because they may belong to a different parent population and exhibit a different evolution 
from the rest of the sources. We also exclude the two GPS sources 0742+103 and 2021+614. In 
our analysis we will use two samples: one containing known FR II sources (hereafter, the "known 
FR II sample" ) and a second sample comprising both known FR II and uncertain FR class sources 
(hereafter, the "full sample"). The "known FR II sample" contains 103 sources (91 quasars, 3 FR II 
gala xies, 1 BL Lac and 8 sources of the intermediate BL Lac/High Polarization Quasar (HPQ) class; 
e.g.. IVeron-Cetty fc Veronll2000l ). One of the sources in this sample lacks redshift information. The 



"full sample" contains 127 sources (91 quasars, 3 FR II galaxies, 10 BL Lacs, 17 BL Lac/HPQs 
and 6 sources without optical counterparts). In this sample 12 sources lack redshift information. 

In Figure [T] we sh ow the luminosity-redshift distribution of sample based on the data of 
Lister &: Homanl (120051 ) . as well as the flux density cutoffs corresponding to the northern and 
southern sky regions. The smallest and largest observed luminosities in our sample are L^ n ~ 
1.19 x 10 25 WHz" 1 and L^ x « 1.03 x 10 29 WHz" 1 and the redshifts range from z min = 0.0491 to 

^max = 3.408. 
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2.2. Missing Redshifts 

Despite considerable observational effort, the redshift information on the MO J AVE sample is 
incomplete, because of the featureless optical spectra of several blazars, and weak/obscured optical 
counterparts. We address this problem by building a pool of redshifts from sources which have 
known redshifts and flux densities within 0.15 Jy of the source with the unknown redshift. We then 
randomly select a redshift from this pool to be used as the redshift for that source. Alternatively, 
one could randomly select redshifts from the entire pool of 102 sources for known FR II sample (115 
for the full sample), however, we chose the former method because of the large range of luminosity 
and redshift spanned by the sample. Because there is only one missing redshift in the known 
FR II sample, the number of possible redshift combinations is only 13 for the known FR II sample, 
compared to 4.3 x 10 13 for the full sample. In the discussion that follows, we use 13 (for the known 
FR II sample) and 1000 (for the full sample) realizations of the randomized redshifts to determine 
the statistical errors on our best fit model parameters arising from missing redshift information. 



3. Method 

3.1. Parameterized Luminosity Function 

The differential luminosity function of a population of objects is defined as the number of 
objects per unit co-moving volume per unit luminosity interval, i.e., 

where N is the number of objects of luminosity C found in the co-moving volume V at red- 
shift z. Studies of flu x-limited AGN samples using the < V/V max > test, including MO J AVE 



(j Arshakian et al.l 120061 ) . indicate that the RLF generally evolves with redshift. Without losing 



generality, we can write the RLF as 

4>{c,z) = M£)fcv(c,z) (2) 

where <f>o(£) is the local [z ~ 0) RLF and / e v(A z) is the evolution function. 

For the intrinsic RLF, we adopt a parameterization in which the local RLF is a simple power 
law of the form 

no / C\ a 

m = \ t\tj ' (3) 

k 0, elsewhere, 
where L* is an arbitrary constant with units of luminosity and no is a normalization constant. In 
this paper we will use L* = 10 27 WHz . 

Traditionally, the evolution (in the simplest cases taken to be luminosity-independent) has 
been parameterized in two popular forms: a power-law evolution of the form (1 + z) k , or an 



-7- 



exponential evolu tion of the form exp[/cr(z)] where t(z) is the look-back time. Other studies (e.g., 
Willott et al.lll998l ) have used 1- or 2-tailed Gaussian redshift dependencies. We were not able to 
successfully fit the MOJAVE data using these parameterizations. In particular, in several cases 
such parameterizations predicted a large spike in the number of low-redshift sources, which is not 
the case for the MOJAVE sample. Instead, we found that a good fit to the data could be obtained 
using the following luminosity-independent density evolution function: 



fev{£,z) = f D (z) 



z m exp 



2() 



a 



(4) 



where m, zq and a are free parameters of the model. Note that this function does not reduce to 
Id(z) = 1 at z = 0; we therefore assume that the model evolution function is valid for a range of 
redshifts z\ < z < Z2- Combining equations ([2]) through our intrinsic model RLF becomes 



z exp 



a 



which is valid over the domain 



C\ < C < £2 and z\ < z < z%. 



(5) 



(6) 



Because the luminous jet material is moving with a speed comparable to c (bulk Lorentz factor 
7 >> 1), its observed monochromatic luminosity will be boosted as 



L = 5 P C, 



(7) 



where C is the luminosity in the rest frame, p = 2 — a ra( j for continuous jet emission, a ra< j is the 
spectral index, and 5 is the kinematic Doppler factor defined as 

-1 



6 



7 - Vt 5 



1 cos I 



(8) 



where 7 = (l — (3 2 ) is the Lorentz factor and = v/c is the velocity of the emitting plasma. If 
the viewing angle to the jet lies within the range 0° < 9 < 90° and 71 < 7 < 72, then the possible 
Doppler factors range from 

^min = 1/72 (9) 



to 



72 + A/72 - 1 



(10) 



If the intrinsic luminosity C is Doppler boosted as in equation ([7]), then the distribution of 
the observed luminosities L will be different from the distribution of the intrinsic luminosities. 



Following the approach used by iListerl (|2003l ). we derive the form of the observed RLF of the 
Doppler beamed sources as 

HL, *) = -° ( — f D (z) / P s (5)5-^+V dS, (11) 
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where P(S) is the probability density function for 5. This model function is valid over the domain 



where 



and 



L\ < L < L2 and z\ < z < Z2, 
Li = C A 



L 2 = ^maxA- 



In equation (jlip . the limits of integration 5\(L) and ^(-L) are given by 
<5i(L) = min |^ max , max j(5 min , (L/£ 2 ) 1/p }} 

and 

5 2 (L) = max|5 min ,min|(5 max , (L/£i) 1/p j j . 
with 5 m i n and <5 max given by equations Q and (fTUj) . The probability density function for 5 is 

r72 P 7 (7) 



P*(<5) = <T 



/(*) vV - 1 



d7, 



(12) 

(13) 
(14) 

(15) 
(16) 

(17) 



where By (7) is t he probability density function for 7 and the lower limit of integration is give n in 
equation (A6) of iListerl (|2003l ). According to the previous results of lLister &: Marscherl (|1997l ). we 
adopt a power-law form of -By (7) with index k: 



B 7 ( 7 ) = C 7 A 



(18) 



for 7x < 7 < 72, where C is a normalization constant. 

For computational purposes we express Ps(5) using beta functions (see Appendix El eq. IA9 

1 1 k\ „ / Ilk 



as 



B 1 



/ 2 (<5)'2' 2 



(19) 



3.2. Maximum Likelihood Method 

From equation (jlip it is apparent that the model parameters (a, m, Zq, a) of the Doppler- 
beamed RLF are the same as the parameters of the intrinsic RLF. Therefore, we can find the 
parameters of the intrinsic RLF by fitting the Doppler -beamed RLF t o the observed data. For 
this purpose we use the maximum likelihood method of lMarshall et al.l (11983!) . which attempts to 
minimize the function S = — 2 ln(Likelihood). The integral in 5 (eq. 2 of lMarshall et al.lll983l ) 
should be equal to the sample size N for a good fit. Therefore, we must minimize 



N 



S(a,m,z ,a) = -2^1n + 2N 



(20) 



i=l 
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and normalize <&(L,z) such that 



N 



22 , dV 
dz — 

dz 



+ /r 



dz 



Zl 



dV 

dz 



max{Li ,Ljj 
L 2 



dL$(L,z) 



max{Li,L . (2)} 



dL$(L,z), 



(21) 



where ./V is the sample size, and /q ~ 6.00912/47T and /q ~ 2.08012/47T are the fractional area of 
the sky available to the survey (in this section the "+" superscript refers to the northern sky area 
while the "— " superscript refers to the southern sky area: 0° < 5 < —20°). In equation (|2ip we 
take into account that in the MOJAVE sample we have two different non-overlapping sky areas, 
each with its own flux density limit: S^~- = 1.5 Jy and S~- = 2.0 Jy. The L"t- (z) and L~-(z) in 

*> mill ** liiin niin \ / mill \ / 

the equation (|2ip are the monochromatic luminosity limits corresponding to the flux density limits 
of the survey: 

L±» = 4^ inJ Di(z)(l + z)-^-\ (22) 

wher e Dj[z] is the lum inosity distance. To minimize S(a, m, zq, <j) we use the "amoeba" algorithm 
from lPress et al.l (1992). 



Other parameters of the model, such as the redshift limits (z\ and Z2), luminosity limits (£\ 
and £2), power law exponent k of the Lorentz factor distribution (eq. fTHJ) and its range of possible 
values [71,72] are taken as fixed a-priori, and are not included in the set of optimized parameters. 
Some of these parameters (e.g., £\ and k) are poorly constrained (for reasons explained at the end 
of this section), while others can be estimated from the data directly, as follows. 

The sources in both the FR II only sample as well as in the full sample span a broad range 
of redshifts from z m [ n = 0.0491 to 

^max — 3.408. Later in this section we show the upper limit of 
the intrinsic luminosity of the parent population to be about £2 = 10 29 WHz -1 . The flux limit 
of the MOJAVE survey would allow the detection of luminous sources with £ > 10 28 WHz^ 1 at 
z = 4 if their Doppler factors are 5 > 2.15. The lack of such sources in the MOJAVE sample at 
redshifts z > 3.4 may be of importance in modeling the RLF of the parent population. Therefore, 
by extending the range of redshifts to z = 4, we allow more freedom in the optimization procedure. 
In addition, this will ensure that we do not exclude the statistical possibility of some objects at 
higher redshift. Because of this we slightly extend the redshift range and set 



0.04 and z 2 = 4. 



(23) 



Extending the upper redshift limit to higher values should not have a large effect on the model 
RLF, since the flux cutoff of th e survey will make the available comoving volume very small at high 
z (see, e.g., IWillott et al.lll998l ). In addition, because of the above mentioned lack of sources above 
z ~ 3.4, the optimization process will constrain the evolution function (eq. 2|) to vanish rapidly at 
larger redshifts. 



Previous studies (see, e.g.. lHoman et al.l 120061 : ICohen et al.ll2007l ) have shown that the appar- 
ent speeds of powerful AGN jets are closely related to their bulk flow velocities and VLBI core 
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properties. Since the MO J AVE sample contains power ful AGN with highly co re-dominated radio 
structures (jCooper et al.l 120071 ) and superluminal jets (jKellermann et al.l 12004 ) . one might expect 
the parent population to have 71 3> 1. However, t he pa rent population likely contains sources with 
much lower jet speeds, and indeed, ICohen et al.l (|2007l ) estimates that the jet sp eed in Cygnus A 
(one of the sources in our subsamples) is 0.59 < (3 < 0.68. Other authors (e.g., IWardle fc Aaron 
19971 ) obtain similar estimates for jet speeds in kiloparsec scales outflows {j3 > 0.6). Assuming no 
strong deceleration of the jets in FR II sources, we adopt the value of /3 m i n = 0.6, or 71 = 1.25, for 
the minimum jet speeds in the parent population. We will di scuss the effe c ts of this choice on our 
model LF in subsection 14.11 Using recent observational data, I Cohen et al.l (|2007l ) find that for the 
MOJAVE sample, -y ma 



32. 



Lister Marscherl (|1997I ) find, using Monte Carlo simulations, that 
a power-law exponent of the Lore ntz factor distribut ion in the range —1.5 < k < —1.75 provide a 



reasonable fit to the CJ-F survey ([Taylor et al.l ll996). a comparable radio-loud blazar sample. In 
this paper we consider the following range of possible Lorentz factors and the exponent k: 



71 = 1-25, 



72 



32 and k 



•1.5. 



(24) 



We can estimate the lower and upper limits for the intrinsic luminosity as follows. First, from 
equations ©, {TO} and jMD we obtain 



<5min = I/72 = 0.031 and 5 max = 72 + J y% - 1 ~ 64 

and we can apply the equation ([7|) to the observed luminosity range in our sample: L^ n « 1.19 x 
10 25 WILT 1 and L^ x w 1.03 x 10 29 WHz -1 . For example, in the extreme case where the range of 
intrinsic luminosities (£) is maximized, we have L^ n = CiS^^ and = jC2<5 mm an< ^ we obtain 
Ci fts 3.31 x 10 21 WHz" 1 and C 2 ~ 1.05 X 10 32 WHz" 1 . In realit y, there is a very l ow probability 



of having such extreme values in the MOJAVE sample (see, e.g., ICohen et al.l 120071 ) . To fine tune 
this range, we initially fit the data using the values given above. We used the parameters of the 
resulting fitted RLF to produce a large population of sources via Monte Carlo simulations. We 
examined the intrinsic luminosity distribution of a simulated flux-limited sample to see if many 
sources had intrinsic luminosities near the value of C\. If all sources were well above this value, 
we adjusted C\ upward incrementally until we obtained a tight fit of the simulated distribution of 
the intrinsic luminosities to the initial range used in that particular step. A similar procedure was 
applied for the upper limit £2- In this manner we found that for the known FR II sample 

C\ = 10 22 ' 2 R3 1.58 x 10 22 WHz -1 (25) 

and 

£ 2 = 10 29 ' 1 « 1.26 x 10 29 WHz -1 (26) 

provided a good fit to the simulated intrinsic luminosity histogram. Similarly, for the full sample 
we find the following limits for the intrinsic luminosities 

d = 10 21 - 6 w 3.98 x 10 21 WHz" 1 (27) 
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and 



C 2 = 10 29 - 2 ps 1.58 x 10 29 WHz" 1 . (28) 
For these Monte Carlo sim ulations we have used the 64-bit random number generato r (RNG) of 



Marsaglia &: Tsang I (J2004J) as many commonly used algorithms (e.g., Numerical Recipes lPress et al . 
19921 ) lack the necessary resolution for generating deviates that a needed to span the wide range of 
luminosities found in equations (125j) - (128p . 



Substituting the above range of intrinsic luminosities into equations (|13p and (|14p . we obtain 
a theoretical range for the observed luminosities: L\ ~ 1.55 x 10 19 WHz -1 and L 2 ~ 5.15 x 
10 32 W Hz" 1 for the known FR II sample and L x ps 3.89 x 10 18 W Hz" 1 and L 2 « 6.49 x 10 32 W Hz" 1 
for the full sample. These ranges are much larger than the observed range of luminosities in the 
MOJAVE sample (see L^ n and 1^ above). Because the model RLF is not well determined 
outside the observed luminosity range, we adopt a conservative approach and adopt the following 
validity range for the luminosities of the observed (beamed) RLF: L\ = 10 25 WHz -1 and L 2 = 
1.1 x 10 29 WHz" 1 . 



Using the values of C\ and C 2 from equations (|25j) and (|26j) . we find that the MOJAVE cutoff 
luminosities are too h i gh for us to observe some important features of the RLF. For example, 
from Fig. 3 of iListerl (120031 ) it is evident that we would need to observe below the luminosity 
L 4 p a 6.49 x 10 25 WHz" 1 for the known FR II sample (L 4 w 1.63 x 10 25 WHz" 1 for full sample) 
(see IListerl 120031 . eq. 9) to probe the region of the RLF that is most susceptible to the changes 
in values of the lower luminosity L\ of the parent population and power-law index k. But in 
our sample we have too few sources with L < L4. For these reasons we chose to estimate some 
parameters of the model from the data as described above, and not to include them in the set of 
optimized parameters. 



4. Results 

4.1. Model Parameters 

Using our adopted form of density evolution (eq. H|) and parameters from equations (j23]), 
(I24p . (125D and (|26j) . we minimized the quantity S(a, m, zq, a) for 1000 (182 for the known FR II 
sample) randomizations of missing redshifts as described in § I2.2L For each fitted parameter, we 
took the median of the distribution as the best fit value of the respective parameter. The best fit 
values of the model RLF thus obtained are presented in Tabl e El The error estim ates for model 



parameters have been obtained using the AS = 1 method (see. lLampton et al.lll976l ). (For a given 
parameter we maximized the likelihood function while keeping all other parameters constant. We 
then varied this parameter until a variation AS" = 1 was obtained.) We have also obtained an 
estimation of the error in the parameter due to missing redshift information at the level of la 
from the values of the parameters for which the fractional cumulative distribution function was 
equal to either 0.683 or 0.317. We found that these errors are negligible compared to the errors 
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computed using the AS = 1 method. We calculated the normalization constant no, space density 
p for L > 1.3 x 10 25 , and parent population K using the best fit values for the model parameters 
a, m,ZQ,a so that equation (|20p yielded the sample size N = 103 for the known FR II sample and 
N = 127 for the full sample. The errors on no and K were also calculated using their cumulative 
distribution functions as described above. We have evaluated the goodness -of-fit of our model RLFs 



using the two-dimensional Kolmogorov-Smirnov (K-S) test as described in lPress et al.1 ()1992l ). The 
Kolmogorov-Smirnov probability (Pks) is a p-value that shows the probability of observing a K-S 
test statistic (Dks) as large or larger than observed one and can be used to reject a model if its 
value is too small. We will accept a model if Pks > 0.2. The K-S probabilities for our models 
indicated good fits to the data with Monte Carlo realizations of missing redshifts, ranging from 
0.67 to 0.79 for the known FR II sample (with D K s from 0.067 to 0.079) and from 0.40 to 0.97 for 
the full sample (with D KS from 0.048 to 0.089). 

From Table [2] we can see that while the slope for the intrinsic luminosity distribution of the 
FR II only sample is slightly shallower than the slope of the full sample, the parameters of the 
density evolution functions of the two samples agree to within la, suggesting that the objects of 
uncertain FR class (9 BL Lacs, 9 BL Lac/HPQ sources and 6 sources without optical classification) 
may actually be of the FR II class. 

The median of average space densities (for L > 1.3 X 10 25 WHz -1 ) computed for 13 (1000 for 
the full sample) randomizations of the unknown redshifts is ~ 1580 Gpc~ 3 for the known FR II 
sample and ~ 4390 Gpc -3 for the full sample. For the full sample the space density is quite large. 
We speculate that this may indicate that the RLF has a different slope for lower luminosities, but 
the lack of low luminosity sources in our sample does not allow us to verify this hypothesis. We 
can also explain this increase in the space density for the full sample by an underestimation of the 
lower intrinsic luminosity C\ (see equation (|27p ) due to the large number of sources with missing 
redshifts in the full sample. 

As discussed in subsection 13.21 the valu e of the lower limit of t he Lorentz factors 71 is not very 



well constrained, with some authors (e.g., lArshakian et al.l 120061 ) suggesting higher values (e.g. 



71 = 3) than the one adopted here (71 = 1.25). Therefore, to investigate whether or not the choice 
of a particular value of 71 has a significant influence on the model RLF, we have repeated the 
computations for the FR Il-only sample using 71 = 3. While we have obtained a slightly different 
intrinsic luminosity range {C\ = 10 2L8 WHz -1 and £2 = 10 29 WHz -1 ), the model RLF parameters 
are essentially unchanged: a = -2.55±0.06, m = 1.48±0.14, z = 1.27±0.09, and a = 0.76±0.10. 
This relative independence of the results on a particular choice of the lower limit 71 for the Lorentz 
factors is d ue to the ins ensitivity of the bright end (L > L4) of the beamed RLF to the values of 71 
and 70 (see Lister! 120031 . Figure 5). However, the fact that the luminosity functions in Figure 5 of 
Listed (120031 ) differ strongly at lower luminosities than L4 means that the choice of 71 potentially 
could modify the predicted parent population sizes. Indeed, for the FR II only sample with 71 = 3, 
we obtain K = (1.7 ± 0.1) x 10 10 , which is about three times larger than the parent population 
predicted by the RLF computed with 71 = 1.25, but this could be due to the lower value of C\ 
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obtained when 71 = 3. 

In Figure[2]we present the integral source counts iV(> S) per unit of solid angle for the observed 
data (known FR II sample), and as predicted by our fitted RLF after it is beamed. The la error 
b ars in thi s and subsequent figures are computed according to Poisson statistics using the method 
of behrelsl (Il98fih . 



4.2. Redshift Distribution 

In Figure [3] we plot the binned redshift distribution and the associated la error bars for the 
known FR II sample (with the missing redshifts replaced with the averages of the "redshift pools" 
as described in § 12. 2p . The solid line represents the predicted redshift distribution for our best fit 
model, while the faint gray lines show the distributions for the 13 randomizations of the missing 
redshifts. We can see that while the missing redshift information creates a tangible uncertainty in 
the redshift distribution, we obtain a reasonably good overall fit to the data. 



4.3. Radio Luminosity Function 



We use the method of iPage &: Carreral (|2000i ) to construct the observed luminosity function. 



In this method, we compute the value of the RLF in a bin with a luminosity interval L m [ n and 
L mm and a redshift interval z m \ n and z max as: 

N 

$est = f 7Tr (29) 

est , 2max , L JV dzdL 

and its uncertainty: 

5N 

<*$est = t 377 . (3°) 

est ,, max ,l dV dzdL 

where N is the number of objects in the bin and 5N its uncertainty. L m ; n (z) is the minimum 
luminosity within the bin at which we can still detect an object. In equations (1291) and (1301) we 



have switched the order of integration compared to the original formulation of IPage fc Carrera 



(2000). 



We use these same equations to compute the binned (i.e., averaged over a luminosity-redshift 
bin) model RLF. We believe this is the most robust way to compare the observed and model RLF 
when binning is involved. We use equations (|2ip to compute the effective "number of sources" N 
in the luminosity-redshift bin of interest (i.e., we replace Z\, Z2, L\ and L2 in eq. [21] with z mm , 
•Zminj L m \ n and L max of the luminosity-redshift bin of interest). 

The fitted and observed RLFs for the MOJAVE sample are presented in the Figure [H It is 
apparent that the averaged model RLF provides a good fit to the sample data. At first glance, it 
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would also appear that both the fitted and observed RLFs obey a broken power-law that steepens at 
higher lu minosities. While a b e amed LF ca n have different slopes for different luminosity intervals 



(see e.g., Illrrv k, Shaferl 11984 ; iListerl l2003l ). this is does not explain the observed break in the 



observed (beamed) LF (e.g., because L4 is too low for our sample; see Figure [5|). Based on this 



Arshakian et al.l 120061 claimed that a double power-law intrinsic LF is needed to describe the 



observed (beamed) LF. However, we find that this is an artifact of the binning method. 

In Figure [5] we plot the Doppler beamed differential model RLF for the known FR II sample as 
well as its average over bins of varying sizes. We can see that our differential Doppler beamed RLF 
is in fact very close to a simple power law RLF, with only a slight flattening for L < L4. Because of 
the large bins, the agreement between the binned RLF (Figure 0]) and the differential RLF (Figure 
[5]) is apparent only at high luminosities. This can be improved by using smaller bin sizes, but at 
the expense of larger Poisson errors. For a steep RLF with strong evolution across a bin, large 
bins can create apparent breaks in the observed RLFs when these bins intersect the luminosity 
cutoff of the sample (see eq. I22j) . This is because the averages of RLFs computed over parts of the 
bins above the luminosity cutoff will be very different from the averages computed over whole bins 
(see Figure [6]). A superior way to plot a binned LF would be to use "centers of mass" of the bins 
over which averaging is done instead of their geometrical centers but unfortunately, this is almost 
impossible to accomplish without a priori knowledge of the luminosity function slope and evolution 
parameters. We can see from Figure [6] that "the center of mass" of the bins is different from the 
geometrical center, even for bins lying entirely above the flux cutoff (i.e., bin "A" in Figure [6]). 
The presence of a redshift dependence in the luminosity function will shift the "center of mass" of 
the bins along the redshift axis as well. The "chopped" bins will not be centered around the same 
redshifts as the whole bins, and therefore, the average of the RLF computed over a "chopped" 
bin should actually belong to a RLF computed at a different cosmological epoch. That is, by 
assuming that the value of the RLF at the center of the bin is equal to the average RLF computed 
over a smaller part of the bin (e.g., the hatched area of the bin "C" in Figure [6j which may even 
not contain the center of the bin) one can introduce a large error in the case of strongly evolving 
functions. In Figure [7] we present a zoomed in version of the Figure [5] on which, additionally, we 
plot the "break" luminosities. We can see that the shifting of the "center of mass" towards lower 
luminosities for "full" bins produce a shift of the LF to the right for logL > logLg + (AlogL)/2. 
For log L < log L\ + (A log L)/2, the LF flattens because of the shifting of the "center of mass" for 
the "chopped" bins towards larger luminosities as well as towards smaller redshifts (see Figure [6]) . 
We conclude that the double power-law that we see in the observed RLF is, therefore, an artifact 
created by the effect of flux density cutoff on steep power law of the intrinsic RLF combined with 
a strong evolution of the luminosity function acr oss a bin. This double power law is not a property 



of the intrinsic RLF of the MOJAVE sample as lArshakian et al.l (|2006l ) concluded. In fact, we are 



able to obtain good fit using a simple power-law intrinsic RLF. 

As previously mentioned in § 13.21 the RLF depends most strongly on the lower luminosity cutoff 
C\ and Lorentz factor distribution power-law index k at luminosities smaller than L4. At larger 
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luminosities it appears more like a featureless simple power-law, as illustrated in Figure Indeed, 
from mathematical considerations, the slope of the beamed RLF is expected to be nearly identical 
to the slope of the unbeamed (intrinsic) RLF for luminosities between La ~ 6.49 x 10 25 WHz -1 



and L 8 « 5 x 10 29 WHz" 1 for the known FR II sample and La 
L 8 w 4.8 x 10 29 WHz" 1 fsee lLister|[20oi : lUrrv et al.lll991n . 



1.63 x 10 25 WHz" 1 and 



I 



5. Conclusions 

We have analyzed the redshift and flux density distributions of a complete sample of AGN 
selected on the basis of relativistically beamed 15 GHz radio flux density (MOJAVE) to derive 
the parent luminosity function of bright radio-loud blazars. We carried out our analysis on two 
samples, one consisting of only the 103 known FR II class radio sources in MOJAVE ("the known 
FR II sample") and a "full sample" that added 24 sources of uncertain FR class. 

1. We find that the observed MOJAVE RLF can be well- fit using a Doppler-boosted, single 
power-law intrinsic RLF with slope a = — 2.53±0.06 for the known FR II sample (a = — 2.65± 
0.06 for the full sample), and a density evolution function of the form z m exp — i ( z ~ z ° ) ' 
with parameters m = 1.4 ± 0.1, zq = 1.29 ± 0.09, and a = 0.76 ± 0.09 for the known FR 
sample (m = 1.6 ± 0.1, zq = 1.18 db 0.09, and a = 0.80 ± 0.1 for the full sample). We assumed 
a power-law Lorentz factor distribution with the exponent k = —1.5 and 1.25 < 7 < 32. Our 
model is valid over the range 0.04 < z < 4 in redshift, 10 22 ' 2 WHz -1 < C < 10 29 - 1 WHz -1 in 
intrinsic luminosity for the known FR II sample (10 21,6 WHz -1 < C < 10 29 2 WHz -1 for the 
full sample), and 10 25 WHz -1 < L < 1.1 X 10 29 WHz -1 in observed luminosity. 

2. We find a good agreement between the fitted RLF parameters of the two samples, suggesting 
that the objects of uncertain FR class in the MOJAVE sample (9 BL Lacs, 9 intermediate 
quasar /BL Lac sources and 6 sources without optical classification) may in fact belong to the 
FR II class. 

3. We have shown that the double power-law shape of the observed (i.e., beamed) MOJAVE 
RLF is an artifact due to large changes of the evolving RLF across a bin and its interaction 
with the lower luminosity cutoff of the survey. We find no evidence for a break in the intrinsic 
blazar RLF above 10 25 WHz" 1 at 15 GHz. 
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A. Representation of P$ (5) using incomplete beta-functions 

In order to avoid direct numerical integration, it is convenient to expr ess the P$ (8) thro ugh 
the beta-functions and then use the continued-fraction representation (see iPress et al.l 1 19921 ) for 
fast computation of the integral. 



If 



P 1 { 1 )=C 1 k 



then 

Px (VTl = r)- 2 

If {5) vV - 1 JfV) vV - 1 



P 5 (5) = 5~ 2 / = C5- 2 / ,1 d 1 = C5- 2 G (k, f (5) , 72 ) (Al) 



where we have defined 



G(k, Zl ,z 2 )= r j== — - dj. (A2) 
Jzi V7 - 1 



Making the substitution t = 7 2 , this integral can be rewritten as: 



G(k,z u z 2 ) = i^%^(t-l)-§d 7 (A3) 



/ t — - 1 (t-l a l dt 
i ( „ f 2 k + 1 l \ o ( 2 k + 1 1 



where 



2[ B ^.—.-J-fll^. 



B(z,a,b) = I* t"- 1 (t - l)^ 1 dt 
Jo 



is the incomplete beta function defined for 0<z<l,a>0 and 6 > 0. In our problem z 2 > z\ > 1 
and therefore we represent the beta function through the hyper-geometric function 

B (z, a, b) = a~ l z a 2 F 1 (a, 1 - 6, a + 1; z) (A4) 
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and continue it analytically into the region \z\ > 1 using the formula (see lLandau &: Lifshitall989l . 
eq. (e.6)): 



2 Fi (a,(3,r,z) 

We then obtain: 

B (z, a, b) 



+ 



r( 7 )r(/?-a) ( 1 

— — — -{-z) 2 F X a,a+l-7,a + l-ft- 

r {p) Tin -a) V 
r( 7 )r(a-/3) „ 



(A5) 



i»r( 7 -/3) 



2^1 A/3 + 1 -7,^+1 -«;- ,N > 1. 



crV 2-F1 (a, 1 - 6, a + l;z) 

, a fr(a + i)r(i-o-6) 



a z 



r(i-6)r(i) 



(-z)~ a 2 Fl a,0,a + 6;- + 



r(a + i)r(o + 6-i) 



(-z) 6-1 2-P1 ( 1 - 6, 1 - a - 6, 2 - a - 6; ■ 



T (a) T (a + 6) 



r(a + i) 



(-1) 



_ a r(l-a-fe) 

r(i-b) 



_! T(a + 6-l) /l 



r(a)T(a + 6) V-z 



l-a-6 



2 Fi 1 - a - b, 1 - 6, 2 - a - 6; - 

z 



and therefore 



5 



1 fc 1 
*T~2'2 



1 fc 1 

2 
1 



z?' 2'2 



, ^2m (m = 0, 1,2, . . 



(A6) 



which can be computed using a continued fraction method (see lPress et al.lll992l ). For situations 
when k is close to 2m (m = 0, 1, 2, . . .) we can use the relationship 



B (z, a,b) = B (a, b) - B (1 - z, 6, a) 



to get: 



Finally, 



G(k,z 1 ,z 2 ) 



B 1 



1 1 k 

~z\'2'~2 



B 1 



C 
2^2 



B 1 



1 1 k 

7|'2'"2 



B 1 



1 1 fc 
^'2'~2 



life 



/ 2 (<5)'2' 2 



(A7) 
(A8) 
(A9) 
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Fig. 1. — The luminosity- redshift distribution of the full MO J AVE sample (omitting 4 known FR Is and two 
GPS sources). Only sources with known redshifts are plotted. We use diamonds for known FR II sources and 
open squares for sources of uncertain FR class. The solid line corresponds to the 1.5 Jy flux density cutoff 
for the sources with positive J2000 declinations and the dashed line corresponds to the 2 Jy flux density 
cutoff for the sources with negative declinations, assuming a flat spectral index. 
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Fig. 2.- 



Integral source count of the known FR II sample N(> S) per unit of solid angle. 
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Fig. 3. — Plot of model redshift distributions for each redshift randomization (thin light gray lines), best 
fit model distribution (thick black line) and the observed redshift distribution of the MOJAVE known FR II 
sample (filled circles with error bars corresponding to ler confidence level). 
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Fig. 4. — The observed binned RLF of the MOJAVE known FR II sample (filled circles and lines) and the 
binned model RLF (open squares) for three redshift intervals: 0.04 < z < 0.7 (red, dotted line), 0.7 < z < 1.2 
(green, dot-dashed line), and 1.2 < z < 4 (blue, solid line). The luminosity bin width is AlogL = 0.5. 
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Fig. 5. — Differential model LF Q(L, z) and differential LFs averaged over bins of varying size. The best 
agreement between the differential LF and its binned version is obtained for smaller bin sizes. For large 
bin sizes the binned (or averaged over bins) LF flattens for lower luminosities. The vertical lines show the 
"break" luminosities for each LF: logL b rcak = log-L~ ia (S~ in = 2 Jy,2^) + (AlogL)/2. 
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Fig. 6. — Illustration of distortions induced by binning and flux cutoff on the binned LF. The curved line 
indicates the flux cutoff: L m i n (S m - ln , z). Rectangles represent bins in the log(Luminosity) — rcdshift plane. 
The value of the binned LF is computed as the average of the differential LF over hatched regions of the bins. 
Filled circles represent the position of the center of the bins (the luminosity coordinates of these centers were 
used to plot the binned LF in Figure [4] and in Figure [5]). The asterisks show the position of the "center of 
mass" of the hatched regions for the bins assuming a &(L, z) — L~ 2 5 . Because of the power-law distribution 
of luminosities (i -2 5 ), the center of mass for the bin "A" lies below its geometrical center. If, in addition, 
evolution is present, then the center of mass will shift along the redshift axis as well. 
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Fig. 7. — Differential model LF z) (continuous line) and the differential LF averaged over bins of size 
Az = 0.06 (0.34 < z < 0.4) and AlogL = 0.125 (the binned LF represented by red asterisks in Figure [5]) 
(dashed line and asterisks). Asterisks are placed at the same positions as in Figure The vertical lines show 
the "break" luminosities logLJ + (AlogL)/2 and logi^ + (AlogL)/2 where L* = L+ in {S+ in = 1.5 Jy, z = 
0.4) and L* 2 = L min (S min = 2 Jy, z = 0.4). 
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Table 1. The MO J AVE sample. 
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Table 1 — Continued 
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Table 1 — Continued 
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Table 1 — Continued 



B1950 Name 
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Note. - Column (1): Source B1950 name; Col- 
umn (2): redshift; Column (3): Flux density at 
15 GHz in Jy; Column (4): Fanaroff- Riley class; Col- 
umn (5): Optical class: BL=BL Lac, Q=quasar, G=radio 
galaxy, BL/HPQ=BL Lac / High Polarization Quasar (e.g., 



Veron-Cettv fc Veronll2000h . U=unidentified. 

* Classification based on luminosity at 1.4 GHz when the mor- 
phology was inconclusive. 

* Classification based on luminosity at 1.4 GHz when a kpc- 
scale image was unavailable. 
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Table 2. Best fit model luminosity function parameters. 



Parameter 


FR II only 


All except FR I 


Units 
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-2.53 ±0.06 


-2.65 ±0.06 
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1.6 ±0.1 
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1.18 ±0.09 
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0.8 ±0.1 
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(2.87 ±0.04) x 10" 10 


(2.22 ±0.04) x 10~ 10 


Mpc" 3 
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(1.579 ± 0.008) x 10 3 


(4.39 ± 0.07) x 10 3 


Gpc" 3 


K 


(5.49 ± 0.04) x 10 9 


(1.55 ±0.03) x 10 11 




Ci 


10 22 ' 2 


10 21.6 
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c 2 


1Q 29.1 


1Q 29.2 


WHz^ 1 



Note. — The errors in the parameters a, m, Zq, and a of the model 
LF were computed using the AS = 1 method as described in the 
text. The error estimates on the normalization factor no, space den- 
sity p for L > 1.3 x 10 25 WHz -1 and parent population K were 
calculated using their cumulative distribution functions as described 
in subsection 14. II 



